Author: Amanda Everitt
Began: 8/22/2018
Finished:: 8/23/2018

[Goal]

[Packages Used]

[Results]

Cell Type % of cells Cluster Known Cell Markers
Neuronal 70.52% 0,1,2 Nrgn, Cnih2
Interneuron 11.54% 5,6,8 Dlx1, Dlx2, Arx, Gad2
Astrocyte 4.98% 4 Aldh1l1, Aqp4, Gja1, Mfge8
Microglia 5.31% 3 Aif1, C1qa, C1qb, Ctss
Endothelial 3.55% 7 Cldn5, Ctla2a, Igfbp7, Kdr
Pericyte 2.14% 10,11 Igf2, Col1a1, Col1a2, Dcn
Oligodendrocyte 1.92% 9 Gpr17, Olig2, Pdgfra, Itpr2

1. Load Data

load(paste0(wd, "/outputs/output_01/experiment.aggregate.filtered.Rdata"))
#Normalize -> Log(CPM)
experiment.aggregate <- NormalizeData(object= experiment.aggregate, normalization.method = "LogNormalize", scale.factor = 10000)

2. Exploratory look at data

genes_to_plot = c("Tbr1", "Fezf2","Bcl11b","Bcl11a","Grin2b","Foxp1","Nrgn","Apoe")
plots_list <- list()

for (i in genes_to_plot){
  a <- as.data.frame(Matrix::colSums(experiment.aggregate@data[i, , drop = FALSE]))
  a$sample_type <- "WT"
  a[grep("NULL", rownames(a)),]$sample_type <- "NULL"
  a[grep("HET", rownames(a)),]$sample_type <- "HET"
  colnames(a) <- c("count","sample_type")
  p1<- ggplot(a, aes(x=count, fill=sample_type)) + 
    geom_density(alpha= 0.3) +
    scale_x_continuous(name="Normalized counts per cell") +
    ylim(0,1) + ylab("Density") +
    ggtitle(paste0("Density of ", i," expression"))
  plots_list[[i]] <- p1
}

marrangeGrob(grobs=plots_list, nrow=2, ncol=2)

GenePlot( experiment.aggregate, "nUMI", "nGene", cex.use = 0.5)

3. Identify variable genes

experiment.aggregate <- FindVariableGenes(
  object = experiment.aggregate,
  mean.function = ExpMean,
  dispersion.function = LogVMR,
  x.low.cutoff = 0.125,
  x.high.cutoff = 4,
  y.cutoff = 0.5, do.plot=F)

length(experiment.aggregate@var.genes) 
## [1] 528

4. Linear Dimensional Reduction (sequencing depth, and % mito)

#Scale and Regress data on sequencing depth and % mito
experiment.aggregate <- ScaleData(
  object = experiment.aggregate, #will populate new slot "scale data", does not change "data" slot
  do.scale = TRUE,
  do.center = TRUE, 
  vars.to.regress = c("nUMI", "nGene","percent.mito")) 
## Regressing out: nUMI, nGene, percent.mito
## 
## Time Elapsed:  1.24776516755422 mins
## Scaling data matrix
#Build PCs
experiment.aggregate <- RunPCA(
  object = experiment.aggregate,
  pc.genes = experiment.aggregate@var.genes, #use variable genes
  pcs.compute = 20, 
  maxit = 500)
## [1] "PC1"
##  [1] "Tubb2b"  "Stmn1"   "Stmn2"   "Tubb3"   "Rtn1"    "Fabp7"   "Gap43"  
##  [8] "Map1b"   "Calm2"   "Ncam1"   "Nnat"    "Sox11"   "Snrpn"   "Nrgn"   
## [15] "Gm3764"  "Meis2"   "Mt3"     "Igfbpl1" "Gria2"   "Nrxn3"   "Tubb5"  
## [22] "Sox2"    "Ccnd2"   "Ckb"     "Dbi"     "Atpif1"  "Dpysl2"  "Ptprz1" 
## [29] "Dlx1"    "Dlx2"   
## [1] ""
##  [1] "Cldn5"    "Esam"     "Igfbp7"   "Flt1"     "Vwa1"     "Pglyrp1" 
##  [7] "Itm2a"    "Slc38a5"  "BC028528" "Egfl7"    "Ctla2a"   "AU021092"
## [13] "Slc22a8"  "Ecscr"    "Slc16a1"  "Cd34"     "Kdr"      "Eng"     
## [19] "Pltp"     "Cd93"     "Fn1"      "Ctsh"     "Lmo2"     "Slc40a1" 
## [25] "Gng11"    "Ramp2"    "Car4"     "Slc7a1"   "Klf2"     "Slc7a5"  
## [1] ""
## [1] ""
## [1] "PC2"
##  [1] "Cldn5"    "Esam"     "Vwa1"     "Slc38a5"  "Igfbp7"   "Flt1"    
##  [7] "Pglyrp1"  "Ctla2a"   "Itm2a"    "AU021092" "Slc22a8"  "Kdr"     
## [13] "Slc16a1"  "Mfsd2a"   "Slc7a1"   "Car4"     "Spock2"   "Fn1"     
## [19] "Ly6c1"    "Ramp2"    "Eng"      "Lsr"      "Cd93"     "Slc7a5"  
## [25] "Gng11"    "Cd34"     "Slco1c1"  "Tfrc"     "Ets1"     "Htra3"   
## [1] ""
##  [1] "C1qc"    "C1qa"    "C1qb"    "Csf1r"   "Fcrls"   "Fcer1g"  "Ctss"   
##  [8] "Tyrobp"  "Pld4"    "Laptm5"  "Fcgr3"   "Trem2"   "Cd68"    "Cx3cr1" 
## [15] "Ly86"    "Ltc4s"   "Lgmn"    "Unc93b1" "P2ry12"  "Hexb"    "Lyz2"   
## [22] "Spi1"    "Gpr34"   "Grn"     "Aif1"    "Ctsz"    "Gmfg"    "Ctsd"   
## [29] "Hexa"    "Olfml3" 
## [1] ""
## [1] ""
## [1] "PC3"
##  [1] "Stmn1"         "Igfbpl1"       "Tubb5"         "Dlx1"         
##  [5] "Dlx6os1"       "Dlx2"          "Tubb3"         "Stmn2"        
##  [9] "Sp8"           "Calm2"         "Ptma"          "Atpif1"       
## [13] "Sp9"           "Nrxn3"         "Dlx5"          "Meg3"         
## [17] "Cdca7"         "Sox11"         "Hmgn2"         "Hist1h2ap"    
## [21] "2810417H13Rik" "Map1b"         "Gad2"          "Tiam2"        
## [25] "Hmgb1"         "Arx"           "Top2a"         "Rnd3"         
## [29] "Sox4"          "Actb"         
## [1] ""
##  [1] "Apoe"    "Slc1a3"  "Ptprz1"  "Plpp3"   "Mfge8"   "Ptn"     "Bcan"   
##  [8] "Gja1"    "Aqp4"    "Pla2g7"  "Atp1a2"  "Ednrb"   "Htra1"   "Slc4a4" 
## [15] "Sparcl1" "Fabp7"   "Igfbp2"  "Slc7a10" "Atp1b2"  "Cd9"     "Cspg5"  
## [22] "Tril"    "Gpr37l1" "Lfng"    "Btbd17"  "Cd63"    "Clu"     "Slc6a11"
## [29] "Dbi"     "Vcam1"  
## [1] ""
## [1] ""
## [1] "PC4"
##  [1] "Gap43"   "Nrgn"    "Snrpn"   "Calm1"   "Calm2"   "Ndufa4"  "Uqcrq"  
##  [8] "Stmn2"   "Usmg5"   "Snhg11"  "Pcp4"    "Atpif1"  "Dbi"     "Ache"   
## [15] "Cplx2"   "Reln"    "Opcml"   "Syp"     "Mt3"     "Bhlhe22" "Fabp7"  
## [22] "Tubb2b"  "Spock1"  "Sst"     "Lrfn5"   "Clstn3"  "Scg2"    "Tubb3"  
## [29] "Vgf"     "Stmn1"  
## [1] ""
##  [1] "Igfbpl1"       "Top2a"         "Sp8"           "Zbtb20"       
##  [5] "Dlx2"          "Nfib"          "Hmgb2"         "Dlx6os1"      
##  [9] "Cdca7"         "Hist1h2ap"     "Dlx1"          "Xist"         
## [13] "Nusap1"        "Sp9"           "Meis2"         "2810417H13Rik"
## [17] "Ube2c"         "Ccnb1"         "Hmgn2"         "Sox11"        
## [21] "Tpx2"          "Ccnd2"         "Cenpf"         "Cks2"         
## [25] "Cdca8"         "Birc5"         "Cdca3"         "H2afv"        
## [29] "Cenpa"         "Cdc20"        
## [1] ""
## [1] ""
## [1] "PC5"
##  [1] "Vtn"      "Igf2"     "Col1a2"   "Cald1"    "Col1a1"   "Myl9"    
##  [7] "Ndufa4l2" "Dcn"      "Emp3"     "Rgs5"     "Higd1b"   "Nupr1"   
## [13] "Rarres2"  "Phlda1"   "Fstl1"    "Bgn"      "S100a11"  "Col4a1"  
## [19] "Nbl1"     "Nr2f2"    "Ifitm3"   "Serpinf1" "Zic1"     "Cthrc1"  
## [25] "Lgals1"   "Eva1b"    "Col4a2"   "Cxcl12"   "Ifitm2"   "Gpx8"    
## [1] ""
##  [1] "Cldn5"    "Slco1c1"  "Slc38a5"  "AU021092" "Ctla2a"   "Slc7a5"  
##  [7] "Pglyrp1"  "Slc16a1"  "Mfsd2a"   "Kdr"      "Slc7a1"   "Spock2"  
## [13] "Vwa1"     "Tfrc"     "Ly6c1"    "Lmo2"     "Cd34"     "Slc40a1" 
## [19] "Cd93"     "Egfl7"    "Slc2a1"   "Flt1"     "Nampt"    "Car4"    
## [25] "Htra3"    "Ptprz1"   "Fli1"     "Slc1a4"   "Lsr"      "Slc22a8" 
## [1] ""
## [1] ""
PCAPlot(object = experiment.aggregate,dim.1 = 1, dim.2 = 2 )

VizPCA(object = experiment.aggregate, pcs.use=1:2)

PCHeatmap(
    object = experiment.aggregate, 
    pc.use = 1:12, 
    cells.use = 500, 
    do.balanced = TRUE, 
    label.columns = FALSE,
    use.full = FALSE
)

5. Determine statistically significant PCs

experiment.aggregate <- JackStraw(
    object = experiment.aggregate, 
    num.replicate = 100, 
    num.pc = 20)
#Too slow above so loading here instead of regenerating
load(file=paste0(wd, "/", out_dir, "/experiment.aggregate.norm.RData"))
PCElbowPlot(experiment.aggregate,num.pc = 40)

JackStrawPlot(object = experiment.aggregate, PCs = 1:40, nCol = 5)
## Warning: Removed 16527 rows containing missing values (geom_point).

## An object of class seurat in project Siavash scRNA 
##  14933 genes across 17396 samples.
percentVar <- experiment.aggregate@dr$pca@sdev^2/sum(experiment.aggregate@dr$pca@sdev^2)*100
d<- as.data.frame(percentVar)
d$PC <- as.factor(paste0("PC",rownames(d)))
positions <- d$PC
p<-ggplot(data=d, aes(x=PC, y=percentVar)) +
  geom_bar(stat="identity",fill="steelblue") + 
  theme(axis.text.x=element_text(angle=90, hjust=1), legend.position="none") +
  #geom_text(aes(label=round(percentVar,digits = 1)), vjust=1.6, color="white", size=3.5)+
  xlab("Principle Component") + 
  ylab("Variance Explained (%)") +
  scale_x_discrete(limits = positions)+
  scale_y_continuous(breaks = seq(0, 20, by = 1))+
  geom_hline(yintercept=2, col="red")
p

use.pcs = c(1:11)  #Chose all PCs accounting for more than 2% of variance
experiment.aggregate <- FindClusters(
    object = experiment.aggregate, 
    reduction.type = "pca", 
    dims.use = use.pcs, 
    resolution = seq(0.3,1.2,0.3), 
    print.output = FALSE, 
    force.recalc = TRUE,
    save.SNN = TRUE
)
save("experiment.aggregate", file=paste0(wd, "/", out_dir, "/experiment.aggregate.norm.RData"))
use.pcs = c(1:11)
sapply(grep("^res",colnames(experiment.aggregate@meta.data),value = TRUE),
       function(x) length(unique(experiment.aggregate@meta.data[,x])))
## res.0.3 res.0.6 res.0.9 res.1.2 
##      12      15      18      20

6. Choose which resolution is best

experiment.aggregate <- SetAllIdent(experiment.aggregate, id = "res.0.3") #temporarily chose one
experiment.aggregate <- RunTSNE( object = experiment.aggregate, reduction.use = "pca", dims.use = use.pcs, do.fast = TRUE)
load(file=paste0(out_dir, "/experiment_aggregate.0.3.RData"))
TSNEPlot(object = experiment.aggregate.0.3, group.by="res.0.3", pt.size=0.5, do.label = T)

TSNEPlot(object = experiment.aggregate.0.3, group.by="res.0.6", pt.size=0.5, do.label = T)

#I chose resolution 0.3 here
experiment.aggregate.0.3 <- experiment.aggregate
table(experiment.aggregate.0.3@ident,experiment.aggregate.0.3@meta.data$orig.ident)
experiment.aggregate.0.3 <- BuildClusterTree(experiment.aggregate.0.3, 
                                             pcs.use = use.pcs, do.reorder = F, 
                                             reorder.numeric = F, do.plot=F)
save(experiment.aggregate.0.3, file=paste0(wd, "/", out_dir, "/experiment_aggregate.0.3.RData"))
TSNEPlot(object = experiment.aggregate.0.3, pt.size=0.5, do.label=TRUE)

TSNEPlot(object = experiment.aggregate.0.3, group.by="orig.ident", pt.size=0.5)

PlotClusterTree(experiment.aggregate.0.3)

7. Useful plots

#Silhoutte Lengths-- how close is this clustering?
pcas <- experiment.aggregate.0.3@dr$pca@cell.embeddings[, use.pcs]
my.dist <- dist(pcas)
my.clusters <- as.numeric(as.factor(experiment.aggregate.0.3@meta.data$res.0.3))

clust.col <- hue_pal()(13)
sil <- silhouette(my.clusters, dist = my.dist)
sil.cols <- clust.col[ifelse(sil[,3] > 0, sil[,1], sil[,2])]
sil.cols <- sil.cols[order(-sil[,1], sil[,3])]

plot(sil, main = paste(length(unique(my.clusters)), "clusters"), 
    border=sil.cols, col=sil.cols, do.col.sort=FALSE)

See if the possible outliter cells from scater should be removed (I decide not to)

load(paste0(wd, "/outputs/output_01/possible_outliers.Rdata"))
experiment.aggregate.0.3@meta.data$scater_outlier <- "NO"
experiment.aggregate.0.3@meta.data[rownames(experiment.aggregate.0.3@meta.data) %in%
                                     possible_outliers,]$scater_outlier <- "YES"
table(experiment.aggregate.0.3@meta.data$scater_outlier)
## 
##    NO   YES 
## 16702   694
TSNEPlot(object = experiment.aggregate.0.3, group.by="scater_outlier", pt.size=0.5, do.return = TRUE)

8. Output Marker Genes for each cluster

markers_all.0.3 <- FindAllMarkers(object = experiment.aggregate.0.3, min.pct = 0.25, 
                                  thresh.use = 0.25, only.pos = TRUE, 
                                  test.use = "wilcox")
write.table(markers_all.0.3, file = paste0(wd, "/", out_dir, "/01_marker_genes.csv"), sep=",", col.names = NA)
save("markers_all.0.3", file=paste0(wd, "/", out_dir, "/markers_all.0.3.RData"))
load(paste0(wd, "/", out_dir, "/markers_all.0.3.RData"))
Cell Type % of cells Cluster Known Cell Markers
Neuronal 70.52% 0,1,2 Nrgn, Rorb, Cnih2
Interneuron 11.54% 5,6,8 Dlx1, Dlx2, Arx, Gad2
Astrocyte 4.98% 4 Aldh1l1, Aqp4, Gja1, Mfge8
Microglia 5.31% 3 Aif1, C1qa, C1qb, Ctss
Endothelial 3.55% 7 Cldn5, Ctla2a, Igfbp7, Kdr
Pericyte 2.14% 10,11 Igf2, Col1a1, Col1a2, Dcn
Oligodendrocyte 1.92% 9 Gpr17, Olig2, Pdgfra, Itpr2 
my.col = c("lightblue", "red")
#Neuronal
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Nrgn","Rorb"), cols.use = my.col, nCol = 2, no.legend = F)

#Interneron
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Dlx1","Dlx2","Arx","Gad2"), cols.use = my.col, nCol = 2, no.legend = F)

#Astrocyte
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Aldh1l1","Aqp4","Gja1","Mfge8"), cols.use = my.col, nCol = 2, no.legend = F)

#Microglia
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Aif1","C1qa","C1qb","Ctss"), cols.use = my.col, nCol = 2, no.legend = F)

#Endothelial
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Cldn5","Ctla2a","Igfbp7","Kdr"), cols.use = my.col, nCol = 2, no.legend = F)

#Pericyte
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Igf2","Col1a1","Col1a2","Dcn"), cols.use = my.col, nCol = 2, no.legend = F)

#Oligo
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Gpr17","Olig2","Pdgfra","Itpr2"), cols.use = my.col, nCol = 2, no.legend = F)

9. Relabel clusters with intial DEX results

#TSNE
current.cluster.ids <- c(0:11)
new.cluster.ids <- c("Neuronal", "Neuronal", "Neuronal","Microglia","Astrocyte",
                     "Interneuron","Interneuron","Endothelial","Interneuron",
                     "Oligodendrocytes","Pericyte","Pericyte")
experiment.aggregate.0.3@ident <- plyr::mapvalues(x = experiment.aggregate.0.3@ident, 
                                                  from = current.cluster.ids, to = new.cluster.ids)
p1 <- TSNEPlot(object = experiment.aggregate.0.3, pt.size = 0.1, do.return=TRUE)
p1 <- p1 + ggtitle("tSNE by cell type")
p2 <- TSNEPlot(object = experiment.aggregate.0.3, group.by="orig.ident", pt.size = 0.1, do.return=TRUE, colors.use = c("#F8766D","#7CAE00","skyblue"))
p2 <- p2 + ggtitle("tSNE by genotype")
grid.arrange(p1,p2,ncol=2)

#Cluster Tree
suppressPackageStartupMessages(library(ape))
data.tree <- experiment.aggregate.0.3@cluster.tree[[1]]
data.tree$tip.label <- c("Neuronal", "Neuronal", "Neuronal","Microglia","Astrocyte",
                     "Interneuron","Interneuron","Endothelial","Interneuron",
                     "Oligodendrocytes","Pericyte","Pericyte")
plot.phylo(x = data.tree, font = 1)
i <- c(1,2,3) #Off one bc indexing
nodelabels(pch = 19)
tiplabels(data.tree$tip.label[i], i, adj = 0)

#TSNE  coloring just neuronal clusters
experiment.aggregate.0.3@meta.data$neuronal <- "not.neuronal" 
wt.cells <- grep("WT", rownames(experiment.aggregate.0.3@meta.data), value = T)
null.cells <- grep("NULL", rownames(experiment.aggregate.0.3@meta.data), value = T)
het.cells <- grep("HET", rownames(experiment.aggregate.0.3@meta.data), value = T)

experiment.aggregate.0.3@meta.data[experiment.aggregate.0.3@meta.data$res.0.3 %in% c(0:2) &
                                     rownames(experiment.aggregate.0.3@meta.data) %in% wt.cells,]$neuronal <- "neuronal.WT"
experiment.aggregate.0.3@meta.data[experiment.aggregate.0.3@meta.data$res.0.3 %in% c(0:2) &
                                     rownames(experiment.aggregate.0.3@meta.data) %in% null.cells,]$neuronal <- "neuronal.NULL"
experiment.aggregate.0.3@meta.data[experiment.aggregate.0.3@meta.data$res.0.3 %in% c(0:2) &
                                     rownames(experiment.aggregate.0.3@meta.data) %in% het.cells,]$neuronal <- "neuronal.HET"
TSNEPlot(object = experiment.aggregate.0.3, group.by="neuronal", pt.size=0.5, do.return = TRUE)

10. Output only neuronal cells

neuronal_cells_to_keep <- rownames(experiment.aggregate.0.3@meta.data[experiment.aggregate.0.3@meta.data$res.0.3 %in% c(0,1,2),])
save(list = c("neuronal_cells_to_keep"), file=paste0(wd, "/", out_dir, "/neuronal_cells.RData"))
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin17.7.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS/LAPACK: /usr/local/Cellar/openblas/0.3.6_1/lib/libopenblasp-r0.3.6.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ape_5.3               scales_1.0.0          cluster_2.1.0        
##  [4] dynamicTreeCut_1.63-1 gridExtra_2.3         Seurat_2.3.4         
##  [7] Matrix_1.2-17         cowplot_1.0.0         ggplot2_3.2.1        
## [10] knitr_1.25           
## 
## loaded via a namespace (and not attached):
##   [1] Rtsne_0.15          colorspace_1.4-1    class_7.3-15       
##   [4] modeltools_0.2-22   ggridges_0.5.1      mclust_5.4.5       
##   [7] htmlTable_1.13.2    base64enc_0.1-3     rstudioapi_0.10    
##  [10] proxy_0.4-23        npsurv_0.4-0        flexmix_2.3-15     
##  [13] bit64_0.9-7         codetools_0.2-16    splines_3.6.0      
##  [16] R.methodsS3_1.7.1   lsei_1.2-0          robustbase_0.93-5  
##  [19] zeallot_0.1.0       Formula_1.2-3       jsonlite_1.6       
##  [22] ica_1.0-2           kernlab_0.9-27      png_0.1-7          
##  [25] R.oo_1.22.0         compiler_3.6.0      httr_1.4.1         
##  [28] backports_1.1.5     assertthat_0.2.1    lazyeval_0.2.2     
##  [31] lars_1.2            acepack_1.4.1       htmltools_0.4.0    
##  [34] tools_3.6.0         igraph_1.2.4.1      gtable_0.3.0       
##  [37] glue_1.3.1          RANN_2.6.1          reshape2_1.4.3     
##  [40] dplyr_0.8.3         Rcpp_1.0.2          vctrs_0.2.0        
##  [43] gdata_2.18.0        nlme_3.1-141        iterators_1.0.12   
##  [46] fpc_2.2-3           gbRd_0.4-11         lmtest_0.9-37      
##  [49] xfun_0.10           stringr_1.4.0       lifecycle_0.1.0    
##  [52] irlba_2.3.3         gtools_3.8.1        DEoptimR_1.0-8     
##  [55] MASS_7.3-51.4       zoo_1.8-6           doSNOW_1.0.18      
##  [58] parallel_3.6.0      RColorBrewer_1.1-2  yaml_2.2.0         
##  [61] reticulate_1.13     pbapply_1.4-2       rpart_4.1-15       
##  [64] segmented_1.0-0     latticeExtra_0.6-28 stringi_1.4.3      
##  [67] foreach_1.4.7       checkmate_1.9.4     caTools_1.17.1.2   
##  [70] bibtex_0.4.2        Rdpack_0.11-0       SDMTools_1.1-221.1 
##  [73] rlang_0.4.0         pkgconfig_2.0.3     dtw_1.21-3         
##  [76] prabclus_2.3-1      bitops_1.0-6        evaluate_0.14      
##  [79] lattice_0.20-38     ROCR_1.0-7          purrr_0.3.2        
##  [82] labeling_0.3        htmlwidgets_1.5.1   bit_1.1-14         
##  [85] tidyselect_0.2.5    plyr_1.8.4          magrittr_1.5       
##  [88] R6_2.4.0            snow_0.4-3          gplots_3.0.1.1     
##  [91] Hmisc_4.2-0         pillar_1.4.2        foreign_0.8-72     
##  [94] withr_2.1.2         fitdistrplus_1.0-14 mixtools_1.1.0     
##  [97] survival_2.44-1.1   nnet_7.3-12         tsne_0.1-3         
## [100] tibble_2.1.3        crayon_1.3.4        hdf5r_1.3.0        
## [103] KernSmooth_2.23-15  rmarkdown_1.16      data.table_1.12.4  
## [106] metap_1.1           digest_0.6.21       diptest_0.75-7     
## [109] tidyr_1.0.0         R.utils_2.9.0       stats4_3.6.0       
## [112] munsell_0.5.0